function [xtt_sample] = DurbinKoopman(y,z,a,x0,Sig0,D,F,H,A,Q,R,G)
% Durbin and Koopman Simulation Smoother
%   Based on Durbin and Koopman 2003

k = size(F,1);
T = size(y,1);
n = size(y,2);

if isempty(G), G=eye(k); end
if isempty(D), D = zeros(k,1); end
if isempty(a), a = zeros(T,1); end
if isempty(A), A = zeros(n,1); end
if isempty(z), z = ones(T,1); end

%% Create Fake Observations

x_fake = nan(T,k);
y_fake = nan(T,n);

epsilon_fake = mvnrnd(zeros(T,k),Q);
eta_fake = mvnrnd(zeros(T,n),R);

% Initialise the recursion
x_fake(1,:) = mvnrnd(x0,Sig0);
y_fake(1,:) = x_fake(1,:)*H' + eta_fake(1,:);

for t = 2:T
    x_fake(t,:) = a(t,:) * D' + x_fake(t-1,:)*F' + epsilon_fake(t,:);
    y_fake(t,:) = z(t,:) * A' + x_fake(t,:)*H' + eta_fake(t,:);
end

%% Construct the Artificial Observations y* = y_true - y_fake
y_fake(isnan(y)) = NaN;

%% Kalman Filtering and Smoothing of Fake Observations


%% Slow Version
% 
% [xtt_fake, Sigtt_fake, ~, Sigttminus1_fake,~] = KalmanFilter(y_fake,z,a,x0,Sig0,D,F,H,A,Q,R,G);
% [xtt_s_fake, ~] = KalmanSmoother(xtt_fake,Sigtt_fake,Sigttminus1_fake,F,D,a);
% % 
% [xtt, Sigtt, ~, Sigttminus1,~] = KalmanFilter(y,z,a,x0,Sig0,D,F,H,A,Q,R,G);
% [xtt_s, ~] = KalmanSmoother(xtt,Sigtt,Sigttminus1,F,D,a);
% % 
% xtt_sample = xtt_s - (xtt_s_fake - x_fake);

%% Fast Version

[xtt_fake, Sigtt_fake, ~, Sigttminus1_fake,~] = KalmanFilter(y-y_fake,z,a,zeros(k,1),Sig0,zeros(k,1),F,H,A,Q,R,G);
[xtt_s_fake, ~] = KalmanSmoother(xtt_fake,Sigtt_fake,Sigttminus1_fake,F,zeros(k,1),a);

xtt_sample = xtt_s_fake + x_fake;

end

